//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   Quadruple_Range_2D.hh
 * \author Jeremy Roberts
 * \date   06/19/2011
 * \brief  Quadruple_Range_2D class definition.
 * \note   Copyright (C) 2011 Jeremy Roberts
 */
//---------------------------------------------------------------------------//
// $Rev:: 116                                           $:Rev of last commit
// $Author:: j.alyn.roberts@gmail.com                   $:Author of last commit
// $Date:: 2011-07-05 16:33:08 +0000 (Tue, 05 Jul 2011) $:Date of last commit
//---------------------------------------------------------------------------//

#ifndef QUADRUPLE_RANGE_2D_HH
#define QUADRUPLE_RANGE_2D_HH

#include "angle/Quadrature.hh"

namespace slabtran
{

//===========================================================================//
/*!
 * \class Quadruple_Range_2D
 * \brief 2D Quadruple Range quadrature class.
 *
 * For two-dimensional transport, Abu-Shumays developed several quadrature
 * sets that are in some sense analogs of the DPn quadrature.  The sets he
 * developed are product quadratures, which means that polar and azimuthal
 * angles are independent. To illustrate, consider a function
 * \f$ f(\Omega_x,\Omega_y) \f$. Using a product quadrature, the integral
 * is approximated as
 * \f[
 *    \int_{\Omega} f(\Omega_x,\Omega_y) d\Omega \approx
 *    \sum^{4N_{\phi}}_{i=1} \sum^{N_{\theta}}_{j=1}
 *        w_i w_j f(\sin \theta_j \cos \phi_i, \sin \theta_j \sin \phi_i ) \, ,
 * \f]
 * If a function \f$ f \f$ is a polynomial in \f$\Omega_x \f$
 * and \f$\Omega_y \f$ of degree \f$ L \f$, a quadrature set that can
 * exactly integrate \f$ f \f$ is denoted order \f$ L \f$.
 *
 * Here, we implement the quadruple range quadrature, which Abu-Shumays
 * estimated would be well-suited for ``highly heterogeneous problems with
 * numerous corner singularities''.   The quadrature essentially provides a
 * good approximation for angular fluxes that are represented by distinct
 * polynomials in \f$\Omega_x \f$ and \f$\Omega_y \f$ within each quadrant.
 * We expect such a quadrature to be well-suited for response function
 * generation, where highly anisotropic fluxes will result due to incident
 * flux conditions restricted to single surfaces of a region.
 *
 * The polar angle quadrature is Abu-Shumays' double-range quadrature, labeled
 * \f$ K^d_L \f$, which satisfies
 * \f[
 	 	 \int^{\pi}_{0} d\theta \sin \theta \sin^{j} \! \theta =
 	 	  \int^{1}_{0} x^j \frac{2xdx}{\sqrt{1-x^2}} \, ,
 * \f]
 * where \f$ x = \sin \theta \f$. To represent terms for \f$ j = 0 \ldots L \f$
 * requires \f$ L+1 \f$ equations, sufficient for determining \f$ (L+1)/2 \f$
 * points and \f$ (L+1)/2 \f$ weights.  Note, this quadrature is nothing more that a
 * Gaussian quadrature with weight \f$ w(x)=2x/\sqrt{1-x^2} \f$ over the range
 * \f$ 0 \leq x \leq 1 \f$.  Note furthermore that \f$ L = 2N_{\theta}-1 \f$.
 *
 * The azimuthal angle quadrature is Abu-Shumays' symmetric quadruple-range
 * quadrature, labeled \f$ I^{q0}_{L} \f$.  This quadrature satisfies
 * \f[
        \int^{\pi/2}_{0} d\phi \cos^n \! \phi \sin^m \! \phi
         \approx \sum^{N_{\phi}}_{i=1} w_i \cos^n \! \phi_i \sin^m \! \phi_i \, .
 * \f]
 * subject to the symmetry constraints
 * \f[
        \phi_{N_{\phi}+1-i} = \frac{\pi}{2} - \phi_i \, , \, \, \, \, \,
         i = 1 \ldots \frac{N_{\phi}+1}{2} \,
 * \f]
 * and
 * \f[
       w_{N_{\phi}+1-i} = w_i \, , \, \, \, \, \, i = 1 \ldots \frac{N_{\phi}+1}{2} \, .
 * \f]
 * Given these constraints, the moment equations are not all independent, and it
 * suffices to consider only integrands of form
 * \f[
       \sin^{4k} \!  \phi\, , \,\, \sin^{2k+1} \! \phi \, , \, \,
        \sin^{4k+1} \! \phi \cos  \phi \, , \,\,\,\,\, k = 0,1\ldots \, .
 * \f]
 * T
 * Suppose we select the number of points \f$ N_{\phi}\f$, which requires
 * \f$ 2N_{\phi} \f$ equations.  The symmetry constraints account for half of these equations.
 * The remaining \f$ N_{\phi} \f$ degrees of freedom allow for exact integration of
 * polynomials of degree \f$ L = N_{\phi}-1 \f$.  Thus, for a given \f$ L \f$, we require twice
 * as many azimuthal points as polar points.  When each set has the same \f$ L \f$, they
 * are said to be compatible, i.e. they have the same order of accuracy.  Moreover,
 * it gives a decisive way to choose the number of polar angles for a given choice
 * of the azimuthal angle count.
 *
 * For the quadratures as implemented, we use double-precision values as
 * generated in Maple and Mathematica in tables from which the full set
 * is generated.  In the future, we could also hard code the full
 * set (via Maple's codegen feature).
 *
 * The constructor uses the number of azimuthal angles, represented
 * by the sn_order parameter.
 *
 * References:
 *
 * Abu-Shumays, I.K. <em>Nuclear Science and Engineering</em>
 *     <b>64</b>, 299-316 (1977).
 *
 */
/*!
 * \example angle/test/testQuadruple_Range_2D.cc
 *
 * Test of Quadrature.
 */
//===========================================================================//

class Quadruple_Range_2D : public Quadrature<_2D>
{
  private:
    //! Base class typedef.
    typedef Quadrature<_2D> Base;

  public:
    // Constructor.
    Quadruple_Range_2D(int sn_order, double norm);

    // Display the quadrature.
    void display() const;

    void name() const { std::cout << "Quadruple_Range_2D" << std::endl; } ;

  private:

    // return cos of theta value
    double get_theta(int N, int i, int j);

    // return sine of phi value
    double get_phi(int N, int i, int j);

};

} // end namespace slabtran

#endif // QUADRUPLE_RANGE_2D_HH

//---------------------------------------------------------------------------//
//              end of Quadruple_Range_2D.hh
//---------------------------------------------------------------------------//
